 
      subroutine maxwell(nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, precon, tiny, &
         ncells, ijkcell, itdim, nsp, ibp1, jbp1, kbp1, vvolaxis, divetil, cdlt&
         , sdlt, strait, dx, dy, dz, clite, cntr, dt, t, fourpi, itmax, error &
         , periodicx, bx_ini, &
         c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, &
         c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvol, vol, x,y,z, qom, qdnv&
         , qdnc, curlx, curly, curlz, exmu, eymu, ezmu, rx, ry, rz, qx, qy, qz, qxm1&
         , qym1, qzm1, aqxm1, aqym1, aqzm1, aqx, aqy, aqz, paqx, paqy, paqz, &
         a11, a12, a13, a21, a22, a23, a31, a32, a33, ijkvtmp, qxtilde, qytilde&
         , qztilde, aqxtilde, aqytilde, aqztilde, jxtilde, jytilde, jztilde, &
         paqxm1, paqym1, paqzm1, gdivx, gdivy, gdivz, qxm2, qym2, qzm2, aqxm2, &
         aqym2, aqzm2, paqxm2, paqym2, paqzm2, qxm3, qym3, qzm3, aqxm3, aqym3, &
         aqzm3, paqxm3, paqym3, paqzm3, bxv, byv, bzv, ex0, ey0, ez0, ex, ey, &
         ez, iter,ex_ext,ey_ext,ez_ext,shock) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
      use boundary
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtx 
      integer  :: nvtxkm 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer  :: ncells 
      integer  :: itdim 
      integer  :: nsp 
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: kbp1 
      integer , intent(in) :: itmax 
      integer , intent(out) :: iter 
      real(double)  :: tiny 
      real(double)  :: cdlt 
      real(double)  :: sdlt 
      real(double)  :: strait 
      real(double)  :: dx,dy,dz 
      real(double)  :: clite 
      real(double)  :: cntr 
      real(double)  :: dt 
      real(double)  :: t
      real(double)  :: fourpi 
      real(double) , intent(in) :: error 
      logical , intent(in) :: precon 
      logical  :: periodicx,ex_ext,ey_ext,ez_ext,bx_ini
      integer  :: ijkvtx(*) 
      integer  :: ijkcell(*) 
      integer  :: ijkvtmp(*) 
      logical :: shock
      real(double)  :: vvolaxis(*) 
      real(double)  :: divetil(*) 
      real(double)  :: c1x(*) 
      real(double)  :: c2x(*) 
      real(double)  :: c3x(*) 
      real(double)  :: c4x(*) 
      real(double)  :: c5x(*) 
      real(double)  :: c6x(*) 
      real(double)  :: c7x(*) 
      real(double)  :: c8x(*) 
      real(double)  :: c1y(*) 
      real(double)  :: c2y(*) 
      real(double)  :: c3y(*) 
      real(double)  :: c4y(*) 
      real(double)  :: c5y(*) 
      real(double)  :: c6y(*) 
      real(double)  :: c7y(*) 
      real(double)  :: c8y(*) 
      real(double)  :: c1z(*) 
      real(double)  :: c2z(*) 
      real(double)  :: c3z(*) 
      real(double)  :: c4z(*) 
      real(double)  :: c5z(*) 
      real(double)  :: c6z(*) 
      real(double)  :: c7z(*) 
      real(double)  :: c8z(*) 
      real(double)  :: vvol(*) 
      real(double)  :: vol(*) 
      real(double)  :: x(*)
      real(double)  :: y(*)
      real(double)  :: z(*)
      real(double)  :: qom(*) 
      real(double)  :: qdnv(itdim,*) 
      real(double)  :: qdnc(*) 
      real(double)  :: curlx(*) 
      real(double)  :: curly(*) 
      real(double)  :: curlz(*) 
      real(double)  :: exmu(*) 
      real(double)  :: eymu(*) 
      real(double)  :: ezmu(*) 
      real(double)  :: rx(*) 
      real(double)  :: ry(*) 
      real(double)  :: rz(*) 
      real(double) , intent(inout) :: qx(*) 
      real(double) , intent(inout) :: qy(*) 
      real(double) , intent(inout) :: qz(*) 
      real(double) , intent(inout) :: qxm1(*) 
      real(double) , intent(inout) :: qym1(*) 
      real(double) , intent(inout) :: qzm1(*) 
      real(double) , intent(inout) :: aqxm1(*) 
      real(double) , intent(inout) :: aqym1(*) 
      real(double) , intent(inout) :: aqzm1(*) 
      real(double)  :: aqx(*) 
      real(double)  :: aqy(*) 
      real(double)  :: aqz(*) 
      real(double)  :: paqx(*) 
      real(double)  :: paqy(*) 
      real(double)  :: paqz(*) 
      real(double)  :: a11(*) 
      real(double)  :: a12(*) 
      real(double)  :: a13(*) 
      real(double)  :: a21(*) 
      real(double)  :: a22(*) 
      real(double)  :: a23(*) 
      real(double)  :: a31(*) 
      real(double)  :: a32(*) 
      real(double)  :: a33(*) 
      real(double)  :: qxtilde(*) 
      real(double)  :: qytilde(*) 
      real(double)  :: qztilde(*) 
      real(double)  :: aqxtilde(*) 
      real(double)  :: aqytilde(*) 
      real(double)  :: aqztilde(*) 
      real(double)  :: jxtilde(*) 
      real(double)  :: jytilde(*) 
      real(double)  :: jztilde(*) 
      real(double) , intent(inout) :: paqxm1(*) 
      real(double) , intent(inout) :: paqym1(*) 
      real(double) , intent(inout) :: paqzm1(*) 
      real(double) , intent(inout) :: gdivx(*) 
      real(double) , intent(inout) :: gdivy(*) 
      real(double) , intent(inout) :: gdivz(*) 
      real(double) , intent(inout) :: qxm2(*) 
      real(double) , intent(inout) :: qym2(*) 
      real(double) , intent(inout) :: qzm2(*) 
      real(double) , intent(inout) :: aqxm2(*) 
      real(double) , intent(inout) :: aqym2(*) 
      real(double) , intent(inout) :: aqzm2(*) 
      real(double) , intent(inout) :: paqxm2(*) 
      real(double) , intent(inout) :: paqym2(*) 
      real(double) , intent(inout) :: paqzm2(*) 
      real(double) , intent(inout) :: qxm3(*) 
      real(double) , intent(inout) :: qym3(*) 
      real(double) , intent(inout) :: qzm3(*) 
      real(double) , intent(inout) :: aqxm3(*) 
      real(double) , intent(inout) :: aqym3(*) 
      real(double) , intent(inout) :: aqzm3(*) 
      real(double) , intent(inout) :: paqxm3(*) 
      real(double) , intent(inout) :: paqym3(*) 
      real(double) , intent(inout) :: paqzm3(*) 
      real(double)  :: bxv(*) 
      real(double)  :: byv(*) 
      real(double)  :: bzv(*) 
      real(double)  :: ex0(*) 
      real(double)  :: ey0(*) 
      real(double)  :: ez0(*) 
      real(double)  :: ex(*) 
      real(double)  :: ey(*) 
      real(double)  :: ez(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp2, jbp2, kbp2, n, ijk, k, j, i, k2, ijk2, i2, nvtmp
      real(double) :: jx0, jy0, jz0, lambda, lambm1, jxtildet, jytildet, &
         jztildet, lambm2, lambm3, dnorm, rsum, bottom, botm1, botm2, botm3, &
         dero, dern, derd, wsqx, wsqy, wsqz, wsaqx, wsaqy, wsaqz, rnorm, alfa, &
         dxnorm, xnorm,factor, norm_in, norm_b
      logical :: testgl, precongl 
!-----------------------------------------------
!
!
!
!
      ibp2 = ibp1 + 1 
      jbp2 = jbp1 + 1 
      kbp2 = kbp1 + 1 
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, &
         vvolaxis) 
!
!     calculate A*E0
!
!     ****************************************************************
!
!     CONSTRUCT JACOBI PRECONDITIONER
!
!     ***************************************************************
!
!
      if (precon) then 
!
         call block (nvtx, ijkvtx, nsp, itdim, iwid, jwid, kwid, c1x, c2x, c3x&
            , c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
            c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvolaxis, vol, vvol, clite&
            , dt, qom, qdnv, cntr, bxv, byv, bzv, exmu, eymu, ezmu, a11, a12, &
            a13, a21, a22, a23, a31, a32, a33) 
!
!
         call ludecomp (nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33) 
!
      else 
!
         dnorm = 1. 
!dir$ ivdep
         do n = 1, nvtx 
!
            ijk = ijkvtx(n) 
!
            a11(ijk) = dnorm 
            a12(ijk) = 0. 
            a13(ijk) = 0. 
!cc
            a21(ijk) = 0. 
            a22(ijk) = dnorm 
            a23(ijk) = 0. 
!ccc
            a31(ijk) = 0. 
            a32(ijk) = 0. 
            a33(ijk) = dnorm 
!
         end do 
!
      endif 
!
!     END PRE-CONDITIIONER
!
!     ********************************************************************
!
!     CALCULATE INITIAL RESIDUE
!
!     **********************************************************************
!
	factor=1d0
!	
      call avec (nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, ijkcell&
         , qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt, sdlt&
         , strait, dx, dy, dz, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu, bxv, byv, bzv, cntr, &
         dt, 1., clite, fourpi, ex, ey, ez, rx, ry, rz,factor) 
!
     call norm_inside(ibp1+1,jbp1+1,kbp1+1,rx,ry,rz,periodicx,norm_in)
!
     call bc_sources(t,ibp1+1,jbp1+1,kbp1+1,nsp,dx,dy,dz,dt, &
                       qdnc,bxv,byv,bzv,x,y,z,gdivx,gdivy,gdivz, &
                       vvol,periodicx,ex_ext,ey_ext,ez_ext,bx_ini,shock,factor)
!
     call norm_boundary(ibp1+1,jbp1+1,kbp1+1,gdivx,gdivy,gdivz,norm_b)
!
!     factor=(dx*dy*dz*dt)**2
	factor = (norm_in/(norm_b+1.d-10))
	write(*,*)'maxwell: factor',factor,norm_in,norm_b
	
!
     call avec (nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, ijkcell&
         , qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt, sdlt&
         , strait, dx, dy, dz, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu, bxv, byv, bzv, cntr, &
         dt, 1., clite, fourpi, ex, ey, ez, rx, ry, rz,factor) 
!
!
     call bc_sources(t,ibp1+1,jbp1+1,kbp1+1,nsp,dx,dy,dz,dt, &
                       qdnc,bxv,byv,bzv,x,y,z,gdivx,gdivy,gdivz, &
                       vvol,periodicx,ex_ext,ey_ext,ez_ext,bx_ini,shock,factor)
!
!	gdivx(:itdim)=0d0
!	gdivy(:itdim)=0d0
!	gdivz(:itdim)=0d0


!dir$ ivdep
      do n = 1, nvtx 
!
         ijk = ijkvtx(n) 
!
!
         rx(ijk) = gdivx(ijk) - rx(ijk) 
         ry(ijk) = gdivy(ijk) - ry(ijk) 
         rz(ijk) = gdivz(ijk) - rz(ijk) 
!
      end do 
!
      call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, rx, ry, rz) 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         rx, ry, rz, periodicx) 

!
!
!     END CALCULATION OF RESIDUE
!
!     ****************************************************************
!
      qxtilde(:ibp2*jbp2*kbp2) = ex(:ibp2*jbp2*kbp2) 
      qytilde(:ibp2*jbp2*kbp2) = ey(:ibp2*jbp2*kbp2) 
      qztilde(:ibp2*jbp2*kbp2) = ez(:ibp2*jbp2*kbp2)  
!
!
      rsum = 0.0 
! ijkgmx=0
!dir$ ivdep
      do n = 1, nvtx 
!
         ijk = ijkvtx(n) 
!
         qx(ijk) = 0.0 
         aqx(ijk) = 0.0 
         paqx(ijk) = 0.0 
!
         qy(ijk) = 0.0 
         aqy(ijk) = 0.0 
         paqy(ijk) = 0.0 
!
         qz(ijk) = 0.0 
         aqz(ijk) = 0.0 
         paqz(ijk) = 0.0 
!
         qxm1(ijk) = 0.0 
         aqxm1(ijk) = 0.0 
         paqxm1(ijk) = 0.0 
!
         qym1(ijk) = 0.0 
         aqym1(ijk) = 0.0 
         paqym1(ijk) = 0.0 
!
         qzm1(ijk) = 0.0 
         aqzm1(ijk) = 0.0 
         paqzm1(ijk) = 0.0 
!
         qxm2(ijk) = 0.0 
         aqxm2(ijk) = 0.0 
         paqxm2(ijk) = 0.0 
!
         qym2(ijk) = 0.0 
         aqym2(ijk) = 0.0 
         paqym2(ijk) = 0.0 
!
         qzm2(ijk) = 0.0 
         aqzm2(ijk) = 0.0 
         paqzm2(ijk) = 0.0 
!
         qxm3(ijk) = 0.0 
         aqxm3(ijk) = 0.0 
         paqxm3(ijk) = 0.0 
!
         qym3(ijk) = 0.0 
         aqym3(ijk) = 0.0 
         paqym3(ijk) = 0.0 
!
         qzm3(ijk) = 0.0 
         aqzm3(ijk) = 0.0 
         paqzm3(ijk) = 0.0 
!
!
         rsum = rsum + abs(gdivx(ijk)) + abs(gdivy(ijk)) + abs(gdivz(ijk)) 
!
!
!
      end do 
!
      if (rsum == 0.0) then 
         iter = 0 
         go to 10 
      endif 
!
      bottom = 0.0 
      botm1 = 0.0 
      botm2 = 0.0 
      botm3 = 0.0 
!
!
!     **************  MAIN ITERATION LOOP  **************************
!     begin conjugate gradient interation 
!
      do iter = 1, itmax 
!
!
!     APPLY BLOCK PRE-CONDITIONING
!
!
         call precond (nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33, rx, ry, rz, qxtilde, qytilde, qztilde) 
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, qxtilde, qytilde, qztilde, periodicx) 
!
!
!
!dir$ ivdep
         do n = 1, nvtx 
!
            ijk = ijkvtx(n) 
!
            qxtilde(ijk) = ex(ijk) + qxtilde(ijk) 
            qytilde(ijk) = ey(ijk) + qytilde(ijk) 
            qztilde(ijk) = ez(ijk) + qztilde(ijk) 
!
         end do 
!
         call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, qxtilde, &
            qytilde, qztilde) 
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, qxtilde, qytilde, qztilde, periodicx) 
!
         call avec (nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, &
            ijkcell, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, &
            cdlt, sdlt, strait, dx, dy, dz, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, &
            c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
            c4z, c5z, c6z, c7z, c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu&
            , bxv, byv, bzv, cntr, dt, 1., clite, fourpi, qxtilde, qytilde, &
            qztilde, aqxtilde, aqytilde, aqztilde,factor) 

            do n = 1, nvtx 
!
               ijk = ijkvtx(n) 
!
               aqxtilde(ijk) = rx(ijk) - (gdivx(ijk)-aqxtilde(ijk)) 
               aqytilde(ijk) = ry(ijk) - (gdivy(ijk)-aqytilde(ijk)) 
               aqztilde(ijk) = rz(ijk) - (gdivz(ijk)-aqztilde(ijk)) 
!
               qxtilde(ijk) = qxtilde(ijk) - ex(ijk) 
               qytilde(ijk) = qytilde(ijk) - ey(ijk) 
               qztilde(ijk) = qztilde(ijk) - ez(ijk) 
!
            end do 
!
         lambda = 0.0 
         lambm1 = 0.0 
         lambm2 = 0.0 
         lambm3 = 0.0 
!
!dir$ ivdep
         do n = 1, nvtx 
!
            ijk = ijkvtx(n) 
!
            lambda = lambda + aqxtilde(ijk)*paqx(ijk) + aqytilde(ijk)*paqy(ijk)&
                + aqztilde(ijk)*paqz(ijk) 
!
            lambm1 = lambm1 + aqxtilde(ijk)*paqxm1(ijk) + aqytilde(ijk)*paqym1(&
               ijk) + aqztilde(ijk)*paqzm1(ijk) 
!
            lambm2 = lambm2 + aqxtilde(ijk)*paqxm2(ijk) + aqytilde(ijk)*paqym2(&
               ijk) + aqztilde(ijk)*paqzm2(ijk) 
!
            lambm3 = lambm3 + aqxtilde(ijk)*paqxm3(ijk) + aqytilde(ijk)*paqym3(&
               ijk) + aqztilde(ijk)*paqzm3(ijk) 
!
!
         end do 
!
         lambda = lambda/(bottom + 1.d-30) 
         lambm1 = lambm1/(botm1 + 1.d-30) 
         lambm2 = lambm2/(botm2 + 1.d-30) 
         lambm3 = lambm3/(botm3 + 1.d-30) 
!
!dir$ ivdep
         do n = 1, nvtx 
!
            ijk = ijkvtx(n) 
!
!
            wsqx = qxtilde(ijk) - lambda*qx(ijk) - lambm1*qxm1(ijk) - lambm2*&
               qxm2(ijk) - lambm3*qxm3(ijk) 
            wsqy = qytilde(ijk) - lambda*qy(ijk) - lambm1*qym1(ijk) - lambm2*&
               qym2(ijk) - lambm3*qym3(ijk) 
            wsqz = qztilde(ijk) - lambda*qz(ijk) - lambm1*qzm1(ijk) - lambm2*&
               qzm2(ijk) - lambm3*qzm3(ijk) 
!
            wsaqx = aqxtilde(ijk) - lambda*aqx(ijk) - lambm1*aqxm1(ijk) - &
               lambm2*aqxm2(ijk) - lambm3*aqxm3(ijk) 
            wsaqy = aqytilde(ijk) - lambda*aqy(ijk) - lambm1*aqym1(ijk) - &
               lambm2*aqym2(ijk) - lambm3*aqym3(ijk) 
            wsaqz = aqztilde(ijk) - lambda*aqz(ijk) - lambm1*aqzm1(ijk) - &
               lambm2*aqzm2(ijk) - lambm3*aqzm3(ijk) 
!
            qxm3(ijk) = qxm2(ijk) 
            aqxm3(ijk) = aqxm2(ijk) 
            paqxm3(ijk) = paqxm2(ijk) 
!
            qym3(ijk) = qym2(ijk) 
            aqym3(ijk) = aqym2(ijk) 
            paqym3(ijk) = paqym2(ijk) 
!
            qzm3(ijk) = qzm2(ijk) 
            aqzm3(ijk) = aqzm2(ijk) 
            paqzm3(ijk) = paqzm2(ijk) 
!
            qxm2(ijk) = qxm1(ijk) 
            aqxm2(ijk) = aqxm1(ijk) 
            paqxm2(ijk) = paqxm1(ijk) 
!
            qym2(ijk) = qym1(ijk) 
            aqym2(ijk) = aqym1(ijk) 
            paqym2(ijk) = paqym1(ijk) 
!
            qzm2(ijk) = qzm1(ijk) 
            aqzm2(ijk) = aqzm1(ijk) 
            paqzm2(ijk) = paqzm1(ijk) 
!
            qxm1(ijk) = qx(ijk) 
            aqxm1(ijk) = aqx(ijk) 
            paqxm1(ijk) = paqx(ijk) 
!
            qym1(ijk) = qy(ijk) 
            aqym1(ijk) = aqy(ijk) 
            paqym1(ijk) = paqy(ijk) 
!
            qzm1(ijk) = qz(ijk) 
            aqzm1(ijk) = aqz(ijk) 
            paqzm1(ijk) = paqz(ijk) 
!
            qx(ijk) = wsqx 
            aqx(ijk) = wsaqx 
!
            qy(ijk) = wsqy 
            aqy(ijk) = wsaqy 
!
            qz(ijk) = wsqz 
            aqz(ijk) = wsaqz 
!
         end do 
!
         call precond (nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33, aqx, aqy, aqz, paqx, paqy, paqz) 
!
         call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, paqx, paqy&
            , paqz) 
!
!
!
         rnorm = 0.0 
         alfa = 0.0 
         bottom = 0.0 
         botm1 = 0.0 
         botm2 = 0.0 
         botm3 = 0.0 
!
!
!
!dir$ ivdep
         do n = 1, nvtx 
!
            ijk = ijkvtx(n) 
!
!
            rnorm = rnorm + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk)) 
!
            alfa = alfa + rx(ijk)*paqx(ijk) + ry(ijk)*paqy(ijk) + rz(ijk)*paqz(&
               ijk) 
!
            bottom = bottom + aqx(ijk)*paqx(ijk) + aqy(ijk)*paqy(ijk) + aqz(ijk&
               )*paqz(ijk) 
!
            botm1 = botm1 + aqxm1(ijk)*paqxm1(ijk) + aqym1(ijk)*paqym1(ijk) + &
               aqzm1(ijk)*paqzm1(ijk) 
!
            botm2 = botm2 + aqxm2(ijk)*paqxm2(ijk) + aqym2(ijk)*paqym2(ijk) + &
               aqzm2(ijk)*paqzm2(ijk) 
!
            botm3 = botm3 + aqxm3(ijk)*paqxm3(ijk) + aqym3(ijk)*paqym3(ijk) + &
               aqzm3(ijk)*paqzm3(ijk) 
!
!
         end do 
!
         alfa = alfa/(bottom + 1.d-30) 
!
!dir$ ivdep
            do n = 1, nvtx 
!
               ijk = ijkvtx(n) 
!
               ex(ijk) = ex(ijk) + alfa*qx(ijk) 
               ey(ijk) = ey(ijk) + alfa*qy(ijk) 
               ez(ijk) = ez(ijk) + alfa*qz(ijk) 
!
               rx(ijk) = rx(ijk) - alfa*aqx(ijk) 
               ry(ijk) = ry(ijk) - alfa*aqy(ijk) 
               rz(ijk) = rz(ijk) - alfa*aqz(ijk) 
!
            end do 
!
         dxnorm = 0.0 
         xnorm = 0.0 
!
!dir$ ivdep
            do n = 1, nvtx 
!
               ijk = ijkvtx(n) 
!
               dxnorm = dxnorm + abs(qx(ijk)) + abs(qy(ijk)) + abs(qz(ijk)) 
!
               xnorm = xnorm + abs(ex(ijk)) + abs(ey(ijk)) + abs(ez(ijk)) 
!
            end do 
!
         dxnorm = dxnorm*abs(alfa) 
!
!      write(*,*)'Iter',iter
!      write(*,*)'TEST Conv2','dxnorm=',dxnorm,'xnorm=',xnorm
!      write(*,*)'TEST Conv2','rnorm=',rnorm,'rsum=',rsum
         if (dxnorm<=error*xnorm .and. rnorm<=error*rsum) go to 10 
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, rx, ry, rz, periodicx) 
!
!
!
      end do 
!
      write (*, *) 'maxwell did not converge', rnorm,rsum,dxnorm,xnorm
!      ex(ijkvtx(:nvtx)) = 0. 
!      ey(ijkvtx(:nvtx)) = 0. 
!      ez(ijkvtx(:nvtx)) = 0. 
!	do n=1,nvtx
!        ijk=ijkvtx(n)
!	write(*,*)n,rx(ijk),ry(ijk),rz(ijk)
!	write(*,*)n,gdivx(ijk),gdivy(ijk),gdivz(ijk)
!	enddo
!	stop
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         ex, ey, ez, periodicx) 
!
!
      return  
!     **********************End of Iteration **************************
!
   10 continue 
      write (*, *) 'maxwell converged', iter, rnorm 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         ex, ey, ez, periodicx) 
!
      call axisvec (ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, ex, ey, ez) 
!
!
!     ********************** END MAXWELL **********************
!
      return  
      end subroutine maxwell 


      subroutine avec(nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, &
         ijkcell, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt&
         , sdlt, strait, dx,dy,dz, periodicx,c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x&
         , c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z&
         , c7z, c8z, vol, curlvx, curlvy, curlvz, exmu, eymu, ezmu, bxv, byv, &
         bzv, cntr, dt, transpos, clite, fourpi, ex, ey, ez, hx, hy, hz,factor) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtmp 
      integer  :: nvtxkm 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer  :: ncells 
      integer  :: nsp 
      integer  :: itdim 
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: kbp1 
      real(double)  :: tiny 
      real(double)  :: factor
      real(double)  :: cdlt 
      real(double)  :: sdlt 
      real(double)  :: strait 
      real(double)  :: dx,dy,dz 
      real(double) , intent(in) :: cntr 
      real(double)  :: dt 
      real(double)  :: transpos 
      real(double)  :: clite 
      real(double)  :: fourpi 
      logical  :: periodicx
      integer  :: ijkvtmp(*) 
      integer  :: ijkcell(*) 
      real(double)  :: qom(*) 
      real(double)  :: qdnv(itdim,*) 
      real(double) , intent(in) :: vvol(*) 
      real(double) , intent(in) :: vvolaxis(*) 
      real(double)  :: c1x(*) 
      real(double)  :: c2x(*) 
      real(double)  :: c3x(*) 
      real(double)  :: c4x(*) 
      real(double)  :: c5x(*) 
      real(double)  :: c6x(*) 
      real(double)  :: c7x(*) 
      real(double)  :: c8x(*) 
      real(double)  :: c1y(*) 
      real(double)  :: c2y(*) 
      real(double)  :: c3y(*) 
      real(double)  :: c4y(*) 
      real(double)  :: c5y(*) 
      real(double)  :: c6y(*) 
      real(double)  :: c7y(*) 
      real(double)  :: c8y(*) 
      real(double)  :: c1z(*) 
      real(double)  :: c2z(*) 
      real(double)  :: c3z(*) 
      real(double)  :: c4z(*) 
      real(double)  :: c5z(*) 
      real(double)  :: c6z(*) 
      real(double)  :: c7z(*) 
      real(double)  :: c8z(*) 
      real(double)  :: vol(*) 
      real(double)  :: curlvx(*) 
      real(double)  :: curlvy(*) 
      real(double)  :: curlvz(*) 
      real(double)  :: exmu(*) 
      real(double)  :: eymu(*) 
      real(double)  :: ezmu(*) 
      real(double)  :: bxv(*) 
      real(double)  :: byv(*) 
      real(double)  :: bzv(*) 
      real(double)  :: ex(*) 
      real(double)  :: ey(*) 
      real(double)  :: ez(*) 
      real(double) , intent(out) :: hx(*) 
      real(double) , intent(out) :: hy(*) 
      real(double) , intent(out) :: hz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk, n, ncellstmp
	  integer, dimension(itdim) :: ijkcelltmp
      real(double) :: rvvol, dtsq, dtc, dtcsq 
	  real(double), dimension(itdim) :: divD
	  logical :: cartesian
!-----------------------------------------------
!
      cartesian = .TRUE.
!
!     calculate the curl of E at cell centers and store in curlc
!
!     AVEC calls CULRC, CURLV, and MUDOTE
!
      call list (1, ibp1, 1, jbp1, 1, kbp1, iwid, jwid, kwid, nvtmp, ijkvtmp) 
!
!
!     EX COMPONENT
!
      call gradc (nvtmp, ijkvtmp, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ex, exmu, eymu, ezmu) 
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvx) 
!
!     ey component
!
      call gradc (nvtmp, ijkvtmp, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ey, exmu, eymu, ezmu) 
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvy) 
!
!
!     EZ COMPONENT
!
      call gradc (nvtmp, ijkvtmp, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ez, exmu, eymu, ezmu) 
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvz) 
!
      call list (2, ibp1 + 1, 2, jbp1 + 1, 2, kbp1 + 1, iwid, jwid, kwid, nvtmp&
         , ijkvtmp) 
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         curlvx, curlvy, curlvz, periodicx) 
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, curlvx, curlvy, curlvz) 
!
!     normalize axis divergence
!
      if (.not.cartesian) then 
         do j = 2, jbp1 + 1 
            do i = 2, ibp1 + 1 
!
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + kwid 
               rvvol = 1./vvol(ijk) 
               curlvx(ijk) = curlvx(ijk)*vvolaxis(ijk)*rvvol 
               curlvy(ijk) = curlvy(ijk)*vvolaxis(ijk)*rvvol 
               curlvz(ijk) = curlvz(ijk)*vvolaxis(ijk)*rvvol 
!
!
            end do 
         end do 
      endif 
!
!
!
!     project susceptibility tensor on to E
!
      call mudote (nvtmp, ijkvtmp, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt, &
         transpos, clite, ex, ey, ez, exmu, eymu, ezmu) 
!
      call graddiv (ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1, kbp1, &
         cdlt, sdlt, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, exmu, eymu, ezmu, divD, hx, hy, hz) 
!
      dtsq = dt**2 
      dtc = clite*cntr*dt 
      dtcsq = dtc**2 
!
	hx(:itdim)=0d0
	hy(:itdim)=0d0
	hz(:itdim)=0d0
!
     if (cartesian) then 
!dir$ ivdep
         do n = 1, nvtmp 
!
 
            ijk = ijkvtmp(n) 
!
            hx(ijk) = (-dtcsq*(curlvx(ijk)+0.5*cntr*dtsq*hx(ijk))) + vvol(ijk)*&
               (0.5*cntr*dtsq*exmu(ijk)+ex(ijk)) 
            hy(ijk) = (-dtcsq*(curlvy(ijk)+0.5*cntr*dtsq*hy(ijk))) + vvol(ijk)*&
               (0.5*cntr*dtsq*eymu(ijk)+ey(ijk)) 
            hz(ijk) = (-dtcsq*(curlvz(ijk)+0.5*cntr*dtsq*hz(ijk))) + vvol(ijk)*&
               (0.5*cntr*dtsq*ezmu(ijk)+ez(ijk)) 
!
         end do 
!
      else 
!dir$ ivdep
         do n = 1, nvtmp 
!
 
            ijk = ijkvtmp(n) 
!
            hx(ijk) = (-dtcsq*curlvx(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*exmu(&
               ijk)+ex(ijk)) 
            hy(ijk) = (-dtcsq*curlvy(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*eymu(&
               ijk)+ey(ijk)) 
            hz(ijk) = (-dtcsq*curlvz(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*ezmu(&
               ijk)+ez(ijk)) 
!
         end do 
      endif 
!
	exmu(:itdim)=0.5*cntr*dtsq*exmu(:itdim)+ex(:itdim)
	eymu(:itdim)=0.5*cntr*dtsq*eymu(:itdim)+ey(:itdim)
	ezmu(:itdim)=0.5*cntr*dtsq*ezmu(:itdim)+ez(:itdim)


      call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, vol, exmu, eymu, ezmu, divD) 
!
     call bc_maxwell(ibp1+1,jbp1+1,kbp1+1,nsp, dx,dy,dz,dt, exmu, eymu, ezmu, divD, &
        bxv, byv, bzv, ex, ey, ez, hx, hy, hz,vvol,periodicx,factor)
!
      return  
      end subroutine avec 


      subroutine axisgrad(ibp1, jbp1, iwid, jwid, kwid, gradxf, gradyf, gradzf) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      real(double)  :: gradxf(*) 
      real(double)  :: gradyf(*) 
      real(double)  :: gradzf(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk 
      real(double) :: gradx, grady, gradz 
!-----------------------------------------------
!
!
!     calculate gradient at the axis by averaging gradients in the poloidal angle
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
      return  
!
      end subroutine axisgrad 


      subroutine axisvec(ibp1, jbp1, iwid, jwid, kwid, vvol, vvolaxis, vecx, &
         vecy, vecz) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ibp1 
      integer  :: jbp1 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      real(double)  :: vvol(*) 
      real(double)  :: vvolaxis(*) 
      real(double)  :: vecx(*) 
      real(double)  :: vecy(*) 
      real(double)  :: vecz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk 
      real(double) :: avgx, avgy, avgz, rvvol 
!-----------------------------------------------
!
!     calculate average vector at the axis by volume averaging
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
      return  
!
      end subroutine axisvec 


      subroutine axisvol(ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol&
         , vvolaxis) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ibp1 
      integer , intent(in) :: jbp1 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: nvtx 
      integer , intent(in) :: ijkvtx(*) 
      real(double) , intent(in) :: vvol(*) 
      real(double) , intent(in) :: vol(*) 
      real(double) , intent(out) :: vvolaxis(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, j, i, ijk 
!-----------------------------------------------
!
      vvolaxis(ijkvtx(:nvtx)) = vvol(ijkvtx(:nvtx)) 
!
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
!
      do j = 2, jbp1 + 1 
!
         vvolaxis(iwid+1+(j-1)*jwid+kwid:ibp1*iwid+1+(j-1)*jwid+kwid:iwid) = &
            0.125*(vol(iwid+1+(j-1)*jwid+kwid:ibp1*iwid+1+(j-1)*jwid+kwid:iwid)&
            +vol(1+(j-1)*jwid+kwid:(ibp1-1)*iwid+1+(j-1)*jwid+kwid:iwid)+vol(1+&
            (j-2)*jwid+kwid:(ibp1-1)*iwid+1+(j-2)*jwid+kwid:iwid)+vol(1+iwid+(j&
            -2)*jwid+kwid:ibp1*iwid+1+(j-2)*jwid+kwid:iwid)) 
!
      end do 
!
      return  
      end subroutine axisvol 

      subroutine bc_maxwell(nxp,nyp,nzp,nsp, dx,dy,dz, dt,exmu, eymu, ezmu,  divD, &
         bxv, byv, bzv, ex, ey, ez, hx, hy, hz,vvol,periodicx,dtsq)
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
      use boundary
 
! Master method to control application of boundary conditions at each face of the system in 3D
! The conditions are applied to all first and last active node in each direction
! The boundary conditions is imposed as a residual in each boundary node
! GMRES will converge to zero residual corresponding to satisfaction of boundary condtions

      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nxp,nyp,nzp
      integer  :: nsp 
      logical  :: periodicx
      real(double)  :: dx,dy,dz,dt 
      real(double)  :: qom(nsp) 
      real(double)  :: exmu(nxp, nyp, nzp) 
      real(double)  :: eymu(nxp, nyp, nzp) 
      real(double)  :: ezmu(nxp, nyp, nzp) 
      real(double)  :: divD(nxp, nyp, nzp) 
      real(double)  :: bxv(nxp, nyp, nzp) 
      real(double)  :: byv(nxp, nyp, nzp) 
      real(double)  :: bzv(nxp, nyp, nzp) 
      real(double)  :: ex(nxp, nyp, nzp) 
      real(double)  :: ey(nxp, nyp, nzp) 
      real(double)  :: ez(nxp, nyp, nzp) 
      real(double)  :: vvol(nxp, nyp, nzp) 
      real(double) , intent(out) :: hx(nxp, nyp, nzp) 
      real(double) , intent(out) :: hy(nxp, nyp, nzp) 
      real(double) , intent(out) :: hz(nxp, nyp, nzp) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i,j,k
      real(double) :: dtsq,div
!-----------------------------------------------
!-----------------------------------------------
!
!	Boundaries in X
!
!	start
!	dtsq = dt**2 
	i=2
      do j=2,nyp
	  do k=2,nzp
	  if (periodicx) then
	  ! do nothing
	  else!

	  hx(i,j,k)=(exmu(i+1,j,k)-exmu(i,j,k))/dx*vvol(i,j,k)*dtsq
	  hy(i,j,k)=ey(i,j,k)-0d0
	  hz(i,j,k)=ez(i,j,k)-0d0
	  end if	 
	  enddo
	  enddo
!
!	end
	
	i=nxp
      do j=2,nyp
	  do k=2,nzp
	  if (periodicx) then
	  hx(i,j,k)=hx(2,j,k)
  	  hy(i,j,k)=hy(2,j,k)
  	  hz(i,j,k)=hz(2,j,k)
	  else
	  hx(i,j,k)=(exmu(i,j,k)-exmu(i-1,j,k))/dx*vvol(i,j,k)*dtsq
	  hy(i,j,k)=ey(i,j,k)-0d0
	  hz(i,j,k)=ez(i,j,k)-0d0
	  end if	 
	  enddo
	  enddo
!
!	Boundaries in Y
!	periodicy is assumed true
!
!	start
	
	j=1
	! do nothing
	j=nyp
      do i=2,nxp
	  do k=2,nzp
	  hx(i,j,k)=hx(i,2,k)
  	  hy(i,j,k)=hy(i,2,k)
  	  hz(i,j,k)=hz(i,2,k)
	  enddo
	  enddo

!
!	Boundaries in Z
!
!	start

!     dtsq = (dx*dy*dz*dt)**2/1.d0
!     dtsq=(dt)**2	
	k=2
      do i=2,nxp
	  do j=2,nyp
	  if (periodicz) then
	  ! do nothing
	  else
	  if(bcMz.eq.1) then
	  ! apply neumann
	  hx(i,j,k)=ex(i,j,k+1)-ex(i,j,k)
	  else
	  !apply dirichelet
	  hx(i,j,k)=ex(i,j,k)-0d0
	  end if
	  if(bcMz.eq.1) then
	  !apply neumann
	  hy(i,j,k)=ey(i,j,k+1)-ey(i,j,k)
	  else
	  !apply dirichelet
	  hy(i,j,k)=ey(i,j,k)-0d0
	  end if
	     div=sum(divD(max(i-1,2):min(i,nxp-1),max(j-1,2):min(j,nyp-1),k))
	     div=div/(1d0+min(i,nxp-1)-max(i-1,2))/(1d0+min(j,nyp-1)-max(j-1,2))
!	  hz(i,j,k)=(ezmu(i,j,k+1)-ezmu(i,j,k))/dz*vvol(i,j,k)*dtsq
	  hz(i,j,k)=div*vvol(i,j,k)*dtsq
!	  hz(i,j,k)=ez(i,j,k)
	  end if	 
	  enddo
	  enddo
!
!	end
	
	k=nzp
      do j=2,nyp
	  do i=2,nxp
	  if (periodicz) then
	  hx(i,j,k)=hx(i,j,2)
  	  hy(i,j,k)=hy(i,j,2)
  	  hz(i,j,k)=hz(i,j,2)
	  else
	  if(bcMz.eq.1) then
	  ! apply neumann
	  hx(i,j,k)=ex(i,j,k)-ex(i,j,k-1)
	  else
	  !apply dirichelet
	  hx(i,j,k)=ex(i,j,k)-0d0
	  end if
	  if(bcMz.eq.1) then
	  !apply neumann
	  hy(i,j,k)=ey(i,j,k)-ey(i,j,k-1)
	  else
	  !apply dirichelet
	  hy(i,j,k)=ey(i,j,k)-0d0
	  end if
	     div=sum(divD(max(i-1,2):min(i,nxp-1),max(j-1,2):min(j,nyp-1),k-1))
	     div=div/(1d0+min(i,nxp-1)-max(i-1,2))/(1d0+min(j,nyp-1)-max(j-1,2))
	  hz(i,j,k)=div*vvol(i,j,k)*dtsq
!	  hz(i,j,k)=(ezmu(i,j,k)-ezmu(i,j,k-1))/dz*vvol(i,j,k)*dtsq
!	  hz(i,j,k)=ez(i,j,k)
	  end if	 
	  enddo
	  enddo

      end subroutine bc_maxwell

      subroutine bc_sources(t,nxp,nyp,nzp,nsp, dx,dy,dz,dt, &
                              qdnc, bxv,byv,bzv,x,y,z,hx, hy, hz, vvol, &
                              periodicx,ex_ext,ey_ext,ez_ext,bx_ini,shock,dtsq)
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
      use boundary
 
! Master method to control application of boundary conditions at each face of the system in 3D
! The conditions are applied to all first and last active node in each direction
! The boundary conditions is imposed as a residual in each boundary node
! GMRES will converge to zero residual corresponding to satisfaction of boundary condtions

      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nxp,nyp,nzp
      integer  :: nsp 
      real(double)  :: dx,dy,dz,dt,t,bx_ini,ex_ext,ey_ext,ez_ext
      logical  :: periodicx,shock
      real(double)  :: qdnc(nxp, nyp, nzp) 
      real(double)  :: vvol(nxp, nyp, nzp) 
      real(double)  :: x(nxp, nyp, nzp) 
      real(double)  :: y(nxp, nyp, nzp) 
      real(double)  :: z(nxp, nyp, nzp) 
      real(double)  :: bxv(nxp, nyp, nzp) 
      real(double)  :: byv(nxp, nyp, nzp)
      real(double)  :: bzv(nxp, nyp, nzp)
      real(double)  :: hx(nxp, nyp, nzp) 
      real(double)  :: hy(nxp, nyp, nzp) 
      real(double)  :: hz(nxp, nyp, nzp) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i,j,k,ibar
      integer :: n,navg
      real(double) :: dtsq,ch,div,pi
      real(double) :: a_m,o_m,g_x,v_int,bavg2,vavg,timescale,z_c
!-----------------------------------------------
!-----------------------------------------------
!
timescale=bx_ini
! the standard for the Newton challenge is a_m=2
!a_m=2.d0*1.5d0
a_m=ey_ext
o_m=.05d0*timescale
ibar=nxp-2
pi=acos(-1.d0)

!     write(*,*)'pre-bc',sum(hz(2:nxp,2:nyp,2:nzp)**2)

!
!	Boundaries in X
!
!	start

!	dtsq = dt**2 
	i=2
      do j=2,nyp
	  do k=2,nzp
	  if (periodicx) then
!	do nothing
	  else
	  hx(i,j,k)=0d0*(qdnc(i,j,k))*vvol(i,j,k)*dtsq
	  hy(i,j,k)=0d0
	  hz(i,j,k)=0d0
	  end if	 
	  enddo
	  enddo
!
!	end
	
	i=nxp
      do j=2,nyp
	  do k=2,nzp
	  if (periodicx) then
	  hx(i,j,k)=hx(2,j,k)
  	  hy(i,j,k)=hy(2,j,k)
  	  hz(i,j,k)=hz(2,j,k)
	  else
	  hx(i,j,k)=0d0*(qdnc(i-1,j,k))*vvol(i,j,k)*dtsq
	  hy(i,j,k)=0d0
	  hz(i,j,k)=0d0
	  end if	 
	  enddo
	  enddo
!
!	Boundaries in Y
!	periodicy is assumed true
!
!	start
	
	j=2
!	do nothing
	j=nyp
      do i=2,nxp
	  do k=2,nzp
	  hx(i,j,k)=hx(i,2,k)
  	  hy(i,j,k)=hy(i,2,k)
  	  hz(i,j,k)=hz(i,2,k)
	  enddo
	  enddo

!
!	Boundaries in Z
!
!	start
!
!        dtsq=(dx*dy*dz*dt)**2/1.d0
!        dtsq=(dt)**2
	k=2
      do i=2,nxp
	  do j=2,nyp
          bavg2=bxv(i,j,k)**2+byv(i,j,k)**2+bzv(i,j,k)**2
          g_x=sin(pi*x(i,j,k)/dx/dble(ibar))**2
          v_int=2.d0*a_m*o_m*tanh(o_m*t)/cosh(o_m*t)**2
          if(shock) then
!	special lines for shock problem
          g_x=1.d0
	  v_int=.1d0
!	end special lines forr shock problem
          end if
	  if (periodicz) then
!	do nothing
	  else
	  hx(i,j,k)=ex_ext
	  hy(i,j,k)=v_int*g_x*bavg2/(abs(bxv(i,j,k))+1.d-10)
	     ch=sum(qdnc(max(i-1,2):min(i,nxp-1),max(j-1,2):min(j,nyp-1),k))
	     div=(1d0+min(i,nxp-1)-max(i-1,2))*(1d0+min(j,nyp-1)-max(j-1,2))
	  hz(i,j,k)=ch/div*vvol(i,j,k)*dtsq
	  end if	 
	  enddo
	  enddo
!
!	end
	
	k=nzp
      do j=2,nyp
	  do i=2,nxp
          bavg2=bxv(i,j,k)**2+byv(i,j,k)**2+bzv(i,j,k)**2
          g_x=sin(pi*x(i,j,k)/dx/dble(ibar))**2
          v_int=2.d0*a_m*o_m*tanh(o_m*t)/cosh(o_m*t)**2
          if(shock) then
!	special lines for shock problem
          g_x=1.d0
	  v_int=.1d0
!	end special lines forr shock problem
          end if
	  if (periodicz) then
	  hx(i,j,k)=hx(i,j,2)
  	  hy(i,j,k)=hy(i,j,2)
  	  hz(i,j,k)=hz(i,j,2)
	  else
	  hx(i,j,k)=ex_ext
	  hy(i,j,k)=v_int*g_x*bavg2/(abs(bxv(i,j,k))+1.d-10)
	     ch=sum(qdnc(max(i-1,2):min(i,nxp-1),max(j-1,2):min(j,nyp-1),k-1))
	     div=(1d0+min(i,nxp-1)-max(i-1,2))*(1d0+min(j,nyp-1)-max(j-1,2))
	  hz(i,j,k)=ch/div*vvol(i,j,k)*dtsq
	  end if	 
	  enddo
	  enddo

!      write(*,*)'post-bc',sum(hz(2:nxp,2:nyp,2:nzp)**2)

      end subroutine bc_sources
      
      subroutine norm_inside(nxp,nyp,nzp,rx,ry,rz,periodicx,norm_in)
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
	  use boundary
      implicit none
      integer :: nxp,nyp,nzp,istart,iend,jstart,jend,kstart,kend
      logical :: periodicx
      real(double)  :: rx(nxp, nyp, nzp) 
      real(double)  :: ry(nxp, nyp, nzp) 
      real(double)  :: rz(nxp, nyp, nzp) 
      real(double) :: norm_in   
      
      if(periodicx) then
      istart=2
      iend=nxp-1
      else
      istart=3
      iend=nxp-1
      end if
      
      jstart=2
      jend=nyp-1
      
      if(periodicz) then
      kstart=2
      kend=nzp-1
      else
      kstart=3
      kend=nzp-1
      end if

         
!
!    compute average norm in the interior nodes
! 
        norm_in=sum(rx(istart:iend,jstart:jend,kstart:kend)**2)+ &
        sum(ry(istart:iend,jstart:jend,kstart:kend)**2)+ &
        sum(rz(istart:iend,jstart:jend,kstart:kend)**2)
        norm_in=norm_in/dble((iend-istart+1)*(jend-jstart+1)*(kend-kstart+1))
	write(*,*)'counter',(iend-istart+1)*(jend-jstart+1)*(kend-kstart+1)
     end subroutine norm_inside
     
           subroutine norm_boundary(nxp,nyp,nzp,hhx,hhy,hhz,norm_b)
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
      implicit none
      integer :: nxp,nyp,nzp,i,j,k,ibp1,jbp1,kbp1,counter
      real(double)  :: hhx(nxp, nyp, nzp) 
      real(double)  :: hhy(nxp, nyp, nzp) 
      real(double)  :: hhz(nxp, nyp, nzp) 
      real(double) :: norm_b  
!
!	compute average norm in the boundary nodes
	norm_b = 0.d0
	
	counter=0
	ibp1=nxp-1
	jbp1=nyp-1
	kbp1=nzp-1
	
!	
!	do j = 2, jbp1
!		do k = 3, kbp1
!		        counter=counter+2
!			norm_b = norm_b + hhx(2,j,k)**2 + hhy(2,j,k)**2 + hhz(2,j,k)**2 &
!				        + hhx(ibp1+1,j,k)**2 + hhy(ibp1+1,j,k)**2 + hhz(ibp1+1,j,k)**2
!	        enddo
!	enddo
!	do k = 3, kbp1
!		do i = 3, ibp1
!			counter=counter+2
!			norm_b = norm_b + hhx(i,2,k)**2 + hhy(i,2,k)**2 + hhz(i,2,k)**2 &
!				        + hhx(i,jbp1+1,k)**2 + hhy(i,jbp1+1,k)**2 + hhz(i,jbp1+1,k)**2
!	        enddo
!	enddo
	do i = 3, ibp1
		do j = 2, jbp1
			counter=counter+2
			norm_b = norm_b + hhx(i,j,2)**2 + hhy(i,j,2)**2 + hhz(i,j,2)**2 &
				        + hhx(i,j,kbp1+1)**2 + hhy(i,j,kbp1+1)**2 + hhz(i,j,kbp1+1)**2
	        enddo
	enddo
!
!	do i = 3, ibp1
!		counter=counter+4
!		norm_b = norm_b + hhx(i,2,2)**2 + hhy(i,2,2)**2 + hhz(i,2,2)**2 &
!				+ hhx(i,2,kbp1+1)**2 + hhy(i,2,kbp1+1)**2 + hhz(i,2,kbp1+1)**2 &
!				+ hhx(i,jbp1+1,2)**2 + hhy(i,jbp1+1,2)**2 + hhz(i,jbp1+1,2)**2 &
!				+ hhx(i,jbp1+1,kbp1+1)**2 + hhy(i,jbp1+1,kbp1+1)**2 + hhz(i,jbp1+1,kbp1+1)**2
!	enddo
!	do j = 3, jbp1
!		counter=counter+4
!		norm_b = norm_b + hhx(2,j,2)**2 + hhy(2,j,2)**2 + hhz(2,j,2)**2 &
!				+ hhx(2,j,kbp1+1)**2 + hhy(2,j,kbp1+1)**2 + hhz(2,j,kbp1+1)**2 &
!				+ hhx(ibp1+1,j,2)**2 + hhy(ibp1+1,j,2)**2 + hhz(ibp1+1,j,2)**2 &
!				+ hhx(ibp1+1,j,kbp1+1)**2 + hhy(ibp1+1,j,kbp1+1)**2 + hhz(ibp1+1,j,kbp1+1)**2
!	enddo
!	do k = 3, kbp1
!		counter=counter+4
!		norm_b = norm_b + hhx(2,2,k)**2 + hhy(2,2,k)**2 + hhz(2,2,k)**2 &
!				+ hhx(ibp1+1,2,k)**2 + hhy(ibp1+1,2,k)**2 + hhz(ibp1+1,2,k)**2 &
!				+ hhx(2,jbp1+1,k)**2 + hhy(2,jbp1+1,k)**2 + hhz(2,jbp1+1,k)**2 &
!				+ hhx(ibp1+1,jbp1+1,k)**2 + hhy(ibp1+1,jbp1+1,k)**2 + hhz(ibp1+1,jbp1+1,k)**2
!	enddo
!!	
!	counter=counter+8
!	norm_b = norm_b + hhx(2,2,2)**2 + hhy(2,2,2)**2 + hhz(2,2,2)**2 &
!				+ hhx(ibp1+1,2,2)**2 + hhy(ibp1+1,2,2)**2 + hhz(ibp1+1,2,2)**2 &
!				+ hhx(2,jbp1+1,2)**2 + hhy(2,jbp1+1,2)**2 + hhz(2,jbp1+1,2)**2 &
!				+ hhx(ibp1+1,jbp1+1,2)**2 + hhy(ibp1+1,jbp1+1,2)**2 + hhz(ibp1+1,jbp1+1,2)**2 &
!				+ hhx(2,2,kbp1+1)**2 + hhy(2,2,kbp1+1)**2 + hhz(2,2,kbp1+1)**2 &
!				+ hhx(ibp1+1,2,kbp1+1)**2 + hhy(ibp1+1,2,kbp1+1)**2 + hhz(ibp1+1,2,kbp1+1)**2 &
!				+ hhx(2,jbp1+1,kbp1+1)**2 + hhy(2,jbp1+1,kbp1+1)**2 + hhz(2,jbp1+1,kbp1+1)**2 &
!				+ hhx(ibp1+1,jbp1+1,kbp1+1)**2 + hhy(ibp1+1,jbp1+1,kbp1+1)**2 + hhz(ibp1+1,jbp1+1,kbp1+1)**2	
!!
	norm_b = norm_b/dble(counter)
	write(*,*)'counter',counter
     end subroutine norm_boundary


